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Several methods developed within the Pierre Auger Collaboration for the estimation of 
the muonic component of the Extensive Air Showers observed in the surface Cherenkov 
detectors are described. The results derived from the data show a deficit of muons 
predicted by the current hadronic interactions models at ultra-high energies. 
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1. Introduction 

The Pierre Auger Observatory is a hybrid detector with the aim of studying the 
characteristics of the longitudinal and lateral profiles of extensive air showers (EAS) 
from ultra high energy cosmic rays (UHECR) with a fluorescence detector array 
(ED) and a surface detector array (SD), respectivelylS^. One of the goals is to 
determine the mass composition of the UHECR. An important feature to distinguish 
composition is the muon content of the shower since heavier primaries will generate 
more muons for the same energy. The predicted number of muons depends on the 
hadronic model used. The water-Cherenkov detectors measure only a combination 
of the much larger electromagnetic component and the muonic component reaching 
the ground, but different techniques are used to recover the muonic content. 

We present the results of three different classes of methods: 1) separating the 
muonic and electromagnetic signals based on their different time structures, 2) 
using inclined showers, for which the electromagnetic component has been strongly 
attenuated in crossing the atmosphere and 3) comparing the ground signals with 
the longitudinal profile, for hybrid events seen also with the ED. Data are compared 
with the hadronic interaction models QGSJETII-Odl^ and EPOS-LHC®. 

2. Estimate of the Muonic Signal in Data 

The Cherenkov photons produced by the shower particles in the SD water- 
Cherenkov detectors (WCD) are sampled by three photomultipliers and digitized 
with EADCs in 25 ns bins. The time response of individual particles is the same for 
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muonic and electromagnetic (EM) particles. Nevertheless, there are two features 
that enable us to separate each component: the signal amplitude and the time-of- 
arrival distribution. The energy deposition in WCD for muons (« 240 MeV) is much 
higher than for electrons and photons (« 10 MeV). The number of EM particles 
is, on average, one order of magnitude higher than muons (at relatively low zenith 
angles). The muon signals are peaked and short and the EM signals are smooth 
and elongated. Concerning the time-of-arrival, the muons arrive earlier than the 
EM particles. Sections |2.1| and 2.2 exploit those features. These methods have 


some limitations that gives rise to variances and systematic bias between models 
and primaries (due to muon pile-up, small muon peaks due to corner clipping on 
the WCD, signal fluctuations). The main source of uncertainty comes from the high 
energy photons that can produce a signal similar to muons. Their contribution is 
expected to be less than 10% to 15% for proton- and iron-induced showers in the 
energy and angular ranges of the methods. 


2.1. Smoothing method 

The EM signal (S'em) is much smoother than the muon signal, so smoothing the 
total signal by averaging over Nuns bins, provides an estimation for the electro¬ 
magnetic contribution. The difference to the original trace is assigned to be the 
muonic component S^. The procedure is repeated for a number of iterations Nuer- 
To avoid under- or over-smoothing, the smoothing window, and Nuer are 

optimized with Monte Carlo simulations to minimize the bias and variances. The 
overall bias obtained for muon fractions for different models and primaries is about 
±0.05 and the average resolution is ±0.08^^. 


2.2. Multivariate method 


The basic idea of the method is to combine the characteristics of the EADC signal 
in a equation to estimate the muon fraction /^. The function was chosen to be 

f^ = a + b0 + c/o ,5 ± d9Po + er ( 1 ) 


where 9 is the reconstructed zenith angle of the shower, f is the distance of the 
WCD to the shower axis, / 0.5 is the portion of the signal larger than 0.5 vEivg 


given by / 0.5 = ^ > 0-5}, is 1 if A is true and 0 otherwise, Xj 

is the signal of the j bin in the signal vector trace x = (xi, ...xn)- The Pq function 
is the normalized zero-frequency component of the power spectrum given by, 

-1 


Pn = 




1 


1 ± 


Hx) 


( 2 ) 


where (x) = S/N is the mean of the signal, cr^(x) is the variance of the signal vector 
and (x^) is the second moment. The / 0.5 and Pq are sensitive to large relative 


'Vertical Equivalent Muon 
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fluctuations and short signals. The parameters (a, b, c, d, e) were fitted to Monte 
Carlo simulations. There are other kinds of parametrizations similar to eq but 
this one was chosen with the objective of minimizing the variance and sensitivity 
of the estimator / to composition and models. The overall bias in for the 

different primaries and models is ±0.02 and the average resolution is ±0.08®. 


2.3. Inclined showers 

In very inclined showers, at zenith angles above 60°, the ground signal in the SD is 
mainly from muons, since the other components have been absorbed higher in the 
atmosphere. The ground signal can be fitted to the equation, 

p = NigPfj_^ig{x,y,9,(j)) (3) 

Pn,i 9 is the model prediction!^ normalized to the reference muon density at the 
ground for proton showers of 10^® eV with QGSJetII-03 and FLUKA interaction 
models. Nig is the shower size parameter, that gives the number of muons relative 
to the reference and can be used as an energy estimator for the shower in this zenith 
angle range. We want the muon content, so we compare N^^ = N^/N^^ig with the 
fractions of the total number of muons on the ground at Monte 

Carlo level {N^^ig is the total number of muons for the reference). For several models 
and primaries and get a resolution of 8% and a systematic uncertainty smaller than 
5%. Finally we can look at the parameter which is A^ig corrected with this 
average bias to (seeP). 

3. Application to data 

The three methods described above were applied to the data collected between 01 
January 2004 and 31 December 2012. A deficit in data was found for the number of 
muons predicted from simulations assuming different hadronic interaction models. 
The smoothing and multivariate methods were applied to events with zenith angles 
9 < 60° and energy E S [lO^®'®®, 10^® °^] eV. Only the detectors with a distance 
from the shower axis f S [950,1050] m were used. At A = 10^® eV the resolutions 
on the core position and energy are about 50 m and 12%, and the energy systematic 
uncertainty is 14%. In figure [^a), we have the results as a function of zenith angle 
for these two methods. The average muon signal found in data divided by the muon 
signal in QGSJetll-04, for a proton shower at 10^® eV, at 1000 m and 9 < 60° is: 
1.33 ± 0.02{stat.) ± 0.05(sys.) (multivariate) 

1.31 ± 0.02(stat.) ± 0 . 09 (s 2 /s.) (smoothing) 

For the inclined shower analysis, the data were hybrid events with zenith angles 
62° < 9 < 80° with an adapted version of the FD quality cuts for these zenith 
angled. The results of the parameter as function of energy are plotted in figure 
Bb)- At 10^®eV, the number of muons in the data exceeds the one from proton 
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(iron) simulations with QGSJetII-03 by a factor of 1.8 (1.4). Note that QGSJetII-04 
and EPOS-LHG have 20% and 30% more muons, at this energy, than QGSJETII-03. 



§n 

(a) 



(b) 


Fig. 1. (a) Muon signal measured as a function of zenith angle for the smoothing and multivariate 

method. Values of the muon rescaling for an energy 10^® eV, at 1000 m from shower axis with 
respect to QGSJetII-04. The error bars and the rectangles are the statistical and systematic 
uncertainties. Points are displaced by ±0.5° for visibility, (b) Average value of it,j/(i?pi3/10^®eV) 
as function of the FD energy, with relation to QGSJetII-03. The gray rectangles are the syst. 
uncertainties of The open circles correspond to the if EpD is varied by its syst. uncertainty. 


4. Study of hybrid events 

At the Pierre Auger Observatory, we have collected a large number of hybrid events, 
for which we measured the longitudinal profile (LP) with the FD and the ground 
signal on the SD. For these events it is possible to simulate air showers with matching 
LP, in order to reduce the ground signal fluctuations event-by-event, and compare 
the simulated ground signal with the data. 

In this study the data collected between 01 January and 31 December 2012 were 
also used. The hybrid events were selected using the cuts defined for the energy 
calibration of the SDI^ and only events with 10^®'® < E < 10^® °^ eV were used. 
The showers were simulated using SENECAISI, with FLUKAESIand for QGSJetII-04 
or EPOS-LHG. For each hybrid data event, simulations with same geometry and 
energy were performed until 12 showers matched X^ax within one sigma of the 
real event. From the 12, only the 3 simulated showers that best matched the real 
LP were considered for full simulation and posterior data comparison. This was 
performed for proton and iron primaries and for QGSJetII-04 and EPOS-LHG. In 
figure a) top, one LP for a data event is plotted with a corresponding simulated 
LP for proton and iron. On the bottom is the ground signal for the same data 
event and simulations. It is clear that the real event has a higher signal than the 
two simulated events chosen to represent its longitudinal profile. 

The signal in the SD at 1000 m, 5'(1000), has different fractions coming from muons 
at different zenith angles. Knowing the evolution the EM and muonic fraction in 
S') 1000) with the angle, from simulation (seelii^), we can obtain a global rescaling 
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factor Re and for the EM and muonic components respectively, where S'(IOOO) 
is then given by equation]^ Since S'(IOOO) is a reconstructed parameter of each 
event, the EM and muon components are rescaled before recalculating the detector 
response for each simulated shower. 

<S'data(1000) = Sresc{RE,Rfi)i,j = .Re5'eM.z,j + R^Rt^R (4) 

where a is the energy scaling of the muonic signal (« 0.89 for both models), i runs 
over each event in data and j labels the primary. 

In figure [^b), the results for Re and are plotted for QGSJetII-04 and EPOS- 
LHC, for mixed and pure proton composition. The mixed composition is the mix 
of p, He, N and Fe which best-fits the observed Xmax distribution. 





0.7 0.8 0.9 1 1.1 1.2 1.3 

Re 


(b) 


Fig. 2. Left: (a) On the top, is the longitudinal profile of one of the hybrid events selected with 
two of its matching simulated showers, for proton and iron. On the bottom is the ground signal 
for the same data event and simulations. Right: (b) Results for Re foi" QGSJetII-04 and 

EPOS-LHC, for mixed and pure proton composition. The ellipses and the gray boxes are the 
statistic and systematic uncertainties, respectively. 


The results obtained are written in table Re are very close to one, while 
shows a rescaling of about 30% to 59%, exposing a deficit of muons in the models. 


Table 1. Re results ±stat. ±syst. uncertainties. 


H. Model 



Re 





Rfj. 



QGSJetII-04 p 

1.09 

± 

0.08 

± 

0.09 

1.59 

± 

0.17 

± 

0.09 

QGSJetII-04 mixed 

1.00 

± 

0.08 

± 

0.11 

1.59 

± 

0.18 

± 

0.11 

EPOS-LHC p 

1.04 

± 

0.08 

± 

0.08 

1.45 

± 

0.16 

± 

0.08 

EPOS-LHC mixed 

1.01 

± 

0.07 

± 

0.08 

1.30 

± 

0.13 

± 

0.09 
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5. Discussion 


In the muon analysis at the Pierre Auger Observatory, a muon deficit has been 
observed in the leading hadronic interaction models with respect to data, consistent 


for different methods. The methods in section 2.1 and 2.2 result in a muon signal 
in data around 1.3 higher than QGSJetII-04 for 10^® eV and 9 < 60°. These values 
are compatible with iron primaries in both models, QGSJetII-04 and EPOS-LHG. 
For inclined showers (with 9 > 62°), in section 2.3 data has around 1.8 more 
muons than QGSJetlT03 at 10^® eV. This corresponds to about 1.5 in relation to 
QGSJetII-04. The inclined results are marginally comparable to the prediction for 
iron showers, but the relative number of muons (Rfj,) seems to increase with energy, 
while for the models it decrease. These discrepancies with models could be related 
to an incorrect energy within the 22% systematic uncertainty of the energy scale or 
problems in the simulation of the hadronic and muonic shower components. 

On the hybrid studies, the muon rescaling in 10 EeV air showers {Ecu = 137 TeV) 
is a factor 1.3 to 1.6 larger than predicted using the leading hadronic interaction 
models tuned to fit the LHG data and lower energy accelerator data. However the 
possibility of having a heavy composition, similar or heavier than iron, is in clear 
contradiction with the measurements of the depth of shower maximum^^. 
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